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A METHOD OF MOSAICING ULTRASONIC VOLUMES FOR 
VISUAL SIMULATION 

FT FT .D AND BACKGROUND OF THE INVENTION 

5 The present invention relates to image processing for medical imaging and, 

more particularly, to a method for registering images, especially ultrasound images, of 
a patient, acquired from several viewpoints. 

The problem addr^ed by the present invention is that of mosaicing a large 
volume from a set of small volumes acquired from a 3D ultrasound device. The 

10 formation of large volumes is necessary for an ultrasound visual simulator, such as the 
one described in US Patent Application No. 08/316,841. Trainees using the simulator 
are able to learn how to identify and diagnose a wide range of medical cases, 
including rare pathologies, by operating a simulated ultrasound device on a 
mannequin, without any need for actual patients. The principle of the simulator is 

15 simple but effective. In an off-line preprocess, a volume buffer is generated from real 
ultrasound images, and then the simulated images are generated by slicing the volume 
buffer on-line. Such images can be generated very rapidly, including postprocessing 
enhancements, and can produce images which are, in most cases, indistinguishable 
from real ultrasound images. 

20 The ultrasonic volume buffer has to be big enough, and lias to represent a large 

portion of the mannequin, to permit the unbounded practice of a real life diagnosis. 
Contemporary ultrasound devices do not provide the capability of obtaining the entire 
volume in a single acquisition. This implies that the volume buffer has to be 
reconstructed from several sub-volumes obtained from different viewpoints. The 

25 registration of mono-modal datasets has been extensively investigated in medical 
applications where atlas data are used. However, ultrasonic datasets are far more 
problematic than other medical-modalities, such as CT or MRI, because the ultrasound 
values are significantly more noisy, blurred, deformed, and have view dependent 
variations which are described below. 

30 The basic technique of a mosaicing operation is to align and register two given 

volumes with a significant overlap into a single volume which smoothly combines the 
information from both. The type of registration technique that can be appropriately 
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applied is directly dependent on the type of variation between the two volumes. Thus, 
to design a registration method for ultrasonic volumes, it is necessary to know the 
type of variation exhibited by ultrasonic volumes. 

The typical size of an ultrasound image generated by common ultrasonic 

5 devices is limited to a maximum width of 12-15 cm. The acquisition of a volume is 
therefore reconstructed from a series of two-dimensional slices. There are two main 
methods of collecting the series of slices: a free-hand collection and a mechanical 
collection. In a free-hand collection location and orientation of the slice is tracked 
by a six-degree-of-freedom device. The slices are stored in the volume and the gaps 

10 between the slices are filled by interpolation. This method is impractical for large 
volumes because of the long acquisition time and the inaccuracy caused by the 
numerous unsampled gaps. 

A better approach is to attach the transducer probe to a mechanical motor that 
sweeps the slice along some type of trajectory (e.g., fan, rotation). In particular, with 

15 a parallel sweep, a series of parallel slices are obtained with no inter-slice gaps. The 
desired image resolution may be obtained by sweeping the slice sufficiently slowly. 
Such parallel dense slices provide small volumes of good quality. A series of such 
volumes needs to be collected and assembled to form a large volume. 

The registration of two volumes requires the detection of the changes between 

20 the two images and the design of a transformation that deforms the images in order to 
remove or reduce the variations between them. The source variations can be classified 
into the following three types. 

1 . Directional Variations 

These variations are due to changes in the view point. They cause a 
25 misalignment that can be simply corrected by a rigid transformation. However, as we 
noted above, the acquisition of the same volume from a different view point causes 
other effects that are not compensated for by spatial transformation. For example, 
shadows are cast with a strong correlation with the probe viewing direction. 

2. Volumetric Variations 

30 These are caused by the characteristics of ultrasonic technology, for example, 

depth gain compensation and gain distortions, and the inherently noisy and blurred 
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ultrasound signal. These effects are difficult to model and to remove. One can 
attempt to reduce them by tuning the acquisition parameters. 
3. Geometric Variations 

Geometric deformations are caused by the movements of the body of the 

5 patient during acquisition. Some movements are forced by the acquisition device, 
because the ultrasound probe must have good contact with the body. The human body 
is soft and not flat, and it is rather difficult to maintain contact without causing the 
patient's muscles to contract. 

Another unavoidable deformation is caused by breathing and other natural 

10 behavior of the patient. Periodic deformation (like that of the heart) can be overcome 
by "gating". In gating, the acquisition is synchronized with the period of the 
deformation, and the slices are acquired in phase with the deformation, using 
equipment similar to E.C.G. that monitors heart activity. 

In most cases, the registration of ultrasound images cannot be based on image 

15 features. In ultrasonic data, it would be very difficult for even a human expert, let 
alone an automatic procedure, to identify useful features, because 3D features can 
appear significantly different in two different images. Therefore, a technique for 
registering ultrasound images must be based directly on image pixel values, and not 
on the interpretation of those values. A reasonable strategy would be to base the 

20 registration technique on some "majority". That is, the success of the technique 
depends on the average success over the whole image, even though part of the image 
may not contain any useful information. Such techniques are classified as automatic 
direct registration techniques. One such technique, for two dimensional images, was 
published by Barber (D. C. Barber (1992), Registration of low resolution medical 

25 images, Phys. Med. Biol., vol. 37 no. 7, pp. 1485-1498, which is incorporated by 
reference for all purposes as if fully set forth herein). Barber's technique is not 
directly suitable for use with ultrasound data sets, however, because of the unique 
characteristics of ultrasound data sets noted above. 

There is thus a widely recognized need for, and it would be highly desirable to 

30 have, a method for automatic direct registration of ultrasound images. 
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SUMMARY OF THE INVENTION 

According to the present invention there is provided a method for warping an 
input target image to match an input reference image, both images being at an original 
level of resolution, both images being composed of a plurality of pixels, the method 

5 comprising the steps of: (a) convolving the input target image with a low pass filter, 
thereby obtaining a filtered target image; (b) convolving the input reference image 
with said low pass filter, thereby obtaining a filtered reference image; (c) subsampling 
said filtered target image, thereby obtaining a subsampled target image; (d; 
subsampling said filtered reference image, thereby obtaining a subsampled reference 

10 image; (e) determining an affine transformation that maximizes a normalized cross- 
correlation between said subsampled target image and said subsampled reference 
image; and (f) scaling said affine transformation to apply to the original level of 
resolution, thereby obtaining a scaled factor of said affine transformation. 

According to the present invention there is provided a method for warping an 

15 input target image to match an input reference image, both images being of a 
dimensionality Z), both images being at an original level of resolution, both images 
being composed of a plurality of pixels, the method comprising the steps of: (a) 
subdividing the input target image into a plurality of target subimages; (b) subdividing 
the input reference image into a plurality of reference subimages, there being a one-to- 

20 one correspondence between said reference subimages and said target subimages; (c) 
for each of said target subimages, determining a local irans formation that warps said 
target subimage to match said corresponding reference subimage. 

According to the present invention there is provided a method for blending a 
three-dimensional target image with a three-dimensional reference image, the target 

25 image including a plurality of points, the method comprising the steps of: (a) 
determining at least one seam plane shared by the target image and the reference 
image; (b) for each of said at least one seam plane: (i) extracting an input target slice 
from the target image along said seam plane, (ii) extracting an input reference slice 
from the reference image along said seam plane, (iii) subdividing the input target slice 

30 into a plurality of target subslices, (iv) subdividing the input reference slice into a 
plurality of reference subslices, there being a one-to-one correspondence between said 
reference subslices and said target subslices, (v) for each of said target subslices, 
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determining a local transformation that warps said target subimage lo match said 
corresponding reference subimage, (vi) recursively performing both said subdividing 
and said determining of said local transformation at a plurality of levels of 
localization, said recursion further including the steps of; (A) for each of said target 
5 subslices, composing said local transformation with a prior transformation, thereby 
obtaining a local composed transformation, and (B) forming a weighted sum of said 
local composed transformations, thereby obtaining a global transformation; (c) 
applying a weighted sulv of said at least one global transformation to at least one of 
the points. 

10 The method of the present invention is built on an affine mapping of a target 

image into a reference image. Because this affine mapping is based on the gradients 
of the images, it is referred to herein as the "gradient technique". The naive gradient 
technique is extended herein by two improvements, referred to herein as 
"multiresolution and localization". The derivations herein are 2D derivations: their 

15 extension to 3D and higher dimensions will be obvious to one skilled in the art. 

The gradient technique assumes that the misalignment of the two images is not 
loo large. That is, the starting position of the two images is "good enough", where the 
exact meaning of "good enough" is discussed below. There are other techniques 
which can register two images with arbitrary partial overlapping. Such techniques 

20 usually perform a coarse alignment, possibly by rigid transformation. The present 
invention includes a partial matching method which applies rigid transformation to 
bring the two images to a satisfactory starting position for the fine and elastic parts of 
the registration. 

The naive gradient technique finds an affine operator T that transforms the 
25 coordinates of the target image into the coordinates of the reference image in a way 
that maximizes the overlap of the two images. This technique is based on a Taylor 
expansion of the exact mapping, under the assumption that there is a majority of target 
pixels which are close enough to the corresponding positions in the reference image. 
The Taylor expansion requires, as its input, gradients of the images These gradients 
30 are approximated numerically. Whenever the displacement is large, however, the 
Taylor expansion is no longer valid. Barber suggested a technique to converge to the 
solution by iterations. The series of transformations are computed, where at each step 
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k the solution for 7* is found over an image warped by the previous transformation 
T k ~ l . The transformations must include some filtering and not cause unnecessary loss 
of information. To minimize applying many transformations, Barber suggested 
composing the transformations. Because Barber uses bilinear transformations, which 

5 are not closed under composition, he approximates his composition by another 
bilinear transformation. 

The solution based on iterations assumes that at each step the majority of 
pixels contribute a displacement tow&rdiMhe correct solution. However, in many 
cases in practice, the initial displacement is too large. The present invention includes 

10 a multiresolution method to better handle large misalignments. The transformation is 
computed on a resolution pyramid and the . results from the low resolution 
transformation are used to guide the computation of finer levels. A resolution 
pyramid consists of the original image and a number of copies of it at lower 
resolutions. At lower resolutions, adjacent pixels and local gradients represent large 

15 distances of the original image. Thus, the low resolution transformation is scaled up 
according to image resolution. Because all the transformations at all levels are affine 
transformations, their composition is also an affine transformation. Thus, the final 
image is a transformation of the original image, and the multiresolution method does 
not cause any information loss. 

20 Once the final transformation is computed it is applied to the target image. 

Then, the transformed image is partitioned into four sub-images and the registration 
process is repeated for each of the sub-images to yield four local transformations 
which refine the registration locally. This process is recursively repeated for each 
sub-image and the registration is further refined. The aggregation of the independent 

25 local affine transformations forms a global elastic transformation with a high degree 
of freedom. 

The above registration method is static in the sense that all the pixels of the 
image contribute to the solution. Many pixels may not contain any useful 
information, or their local environment may not have a matching counterpart in the 
30 other image. This situation often occurs in ultrasound images because of hidden parts 
or noise. Therefore, the present invention includes a rejection mechanism to cull 
away pixels that are estimated to be irrelevant to the solution. Such pairs of pixels are 
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not assigned to the linear system of equations. This mechanism significantly 

improves the registration of partially occluded scenes. 

As noted above, the present invention relies heavily on the Taylor expansion 

and the correctness of the gradients. Thus, the images should be continuous and 
5 differential, which is not necessarily the case in many discrete images. However, the 

registration method need not be affected by sharp local variations. Local variations 

may be softened and blurred by a filter as a preprocess to the gradient calculation. 

The filter also removes part of the noise from the images and smoothes ui x This 

improves the exactness of the Taylor expansion and hence the quality of the results. 
10 However, after the transformation is computed over the filtered images, it is applied to 

the original images during the mosaicing stage. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The invention is herein described, by way of example only, with reference to 
] 5 the accompanying drawings, wherein 

FIG. 1 is a diagram of the data flow of multiresolution registration; 
FIG. 2A is a reference image; 

FIG. 2B is a target image to be warped to match the reference image of FIG. 

2A; 

20 FIG. 2C is the result of warping the target image of FIG. 2B without using 

multiresolution; 

FIG. 2D is the result of warping the target image of FIG. 2B using 
multiresolution; 

FIG. 2E is the image of FIG. 2B after one multiresolution step; 
25 FIG. 2F is the image of FIG. 2B after two multiresolution steps; 

FIG. 3A is a reference image; 

FIG. 3B is a target image to be warped to match the reference image of FIG. 

3A; 

FIG. 3C is the result of warping the target image of FIG. 3B using a global 
30 transformation, including multiresolution; 

FIG. 3D is the result of warping the target image of FIG. 3B using one level of 
localization; 
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FIG. 3E is the result of warping the target image of FIG. 3B, using two levels 
of localization; 

FIG. 4 shows the data flow of the localization process; 

FIG. 5 shows schematically the overlap of a target volume and a reference 
5 volume, including two seam planes; 

FIG. 6 shows the volume mosaicing pipeline. 



DESCRIPTION OF THE PREFERRED EMBODIMENTS 

The present invention is of a method for registration of medical images, 
10 particularly ultrasound images. Specifically, the present invention can be used to 
create an internally consistent numerical database of ultrasound images from a 
collection of disparate, overlapping digital ultrasound subimages. 

The principles and operation of medical image registration according to the 
present invention may be better understood with reference to the drawings and the 
15 accompanying description. 

The Gradient Method 
Given two images, a target image / and a reference image J, the goal is to find 
a transformation T such that L(T(I,J)) is maximal, where L is the normalized cross- 
correlation operator 

%{I(x,y)-I){J(x 9 y)-J) 
20 L(I,J) = ^ ; r (1) 

x,y j,y 

and 7 and J are the average values of the target and reference images, respectively. 
In the development it is assumed that the images are smooth differentiate functions. . 
In the algorithm, a discrete version of the derivatives involved is used. 

The transformation T can be defined as a polynomial which maps the x,y 
25 coordinates of / to the w,v coordinates of J: 

u = a n + a l2 x + a l2 y + a^x 2 + a l5 y 2 + a lb xy + ... 

v = a 2X + a 22 x + a 23 y + a 2A x 2 + a ls y 2 + a l( pcy + ... (2) 
In the preferred embodiments of the present invention, an affine transformation is 
used: 

30 w = <7,| + a ]2 x + a ]3 y 
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but other degrees can be used as well. For example, Barber used a bilinear 
transformation. 

To find the coefficients of the affine transformation, a linear system of 
5 equations is constructed from a large set of pairs of points p, e I and q t e J. Those 
point are implicitly described by their displacement, which can be expressed in terms 
of the transformation. 

Let T(x,y) = (w,v). J warps / onto J, it is assumed that ](xy) = 7(w,v) holds 
for the majority of the image pixels. Define dx = u - x and dy = v - y. Then I(x>y) = 
10 + dx, y + dy). Assuming the displacement vector (dx,dy) is small enough to be 

approximated by the linear term of the Taylor expansion, then 

J{x + dx,y + dy) = J{x 9 y) + dxd x J(x*y) + dyd y J(x,y) (4) 
However, because 

dx - u - a- = a } , + (a n - 1)* + Q\jy 
15 dy = v-y = ai X +a 21 x+(a>n-\)y (5) 

it follows that 

I(x,y) = J(x,y) + (fl„ + (a ]2 - \)x + a n y)d x J(x,y) + (a u + </ 22 x + (a 2 , - 1)^)<V 7 (*, J>) 

( 6 ) 

Assuming that Equation (6) holds for all the pixels in the image, the system of 
20 equations can be written in matrix form. Denote by p- t and q } the values of the target 
and reference pixels, respectively. Let 0 be the vector {</,} and A the vector of 
unknowns: 

where a (m compensates for the amplitude scale of the two images: 

25 flooP/ + •» 

The system of equations in Equation (6) can be written in a matrix form: 



D = 



fa_ fa dqj_ fa fa_ 
Pi dx dx y dx dy dy y dy 



(7) 



then, ideally, 
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Q = DA (8) 
The system of equations is overdetermined, because there are many more pixels than 
unknown coefficients in the vector A. Using a least squares technique, A is obtained 
by a pseudo-inverse method: 
5 A = (&Dy ] D'Q , (9) 

where D 1 is the transpose of matrix D. Because the images are discrete, the partial 
derivatives are approximated by forward and backward differences: 
a r 7(^)*(/(x + l lJ ;)-/(x-l,/;v./2 

d y I(x 9 y) * (l(x,y + \) - I(x,y - 1)) 12 (10) 

10 Multiresolution 

In the following it will be shown how it is possible lo exploit multiresolution 
to improve the registration between two images. Multiresolution improves the 
performance of the method of the present invention for images whose initial 
misalignment is rather large. Because the basic method is based on the Taylor 

15 expansion, the displacement values must be small enough for the gradient values to be 
relevant. At lower resolutions, adjacent pixels and local gradients represent large 
distances of the original image. A displacement computed on a low resolution image 
corresponds to a larger displacement on the highest resolution of the original image. 
These larger displacements may yield transformations that compensate for larger 

20 misalignments. However, those are only approximate transformations, because they 
are based on coarse representations of the original images. 

A resolution pyramid consists of the original image and several copies of it at 
lower resolutions. The multiresolution registration is computed on the resolution 
pyramid, and the results from the low resolution transformation arc used to guide the 

25 computation of the finer levels. The low resolution displacements and the 
corresponding transformations are scaled up according to image resolution. 

Before the multiresolution method is elaborated, it first is shown how the A-th 
factor low resolution images are subsampled. Let I k denote the k-\h factor low 
resolution image of J, a convolution of /with a low-pass filter having a kernel H: 

30 /*(*.J') = S"«2l]>-*/2 I "( fa + , ^ + ^('",/) (ID 

A simple box filter H- l/k 2 yields satisfactory results. 



WO 98/25509 



PCT/US97/20950 



11 



Homogeneous coordinates are used to represent points. A 2D point is 
represented by the triple (xj\w), whose corresponding Cartesian coordinates are 
(xAvjVvv). Accordingly, an affine transformation such as Equation (3) may be 
represented by a matrix: 



5 = 



b u b n 
b 2] b 22 



(12) 



0 0 1 

Given an affine transformation B, its/-th scaled factor i/is defined as: 



B f = 



K b n fxt x 

b 2\ b 22 f X 



(13) 



0 0 1 

and it is easy to verify that if Bp = </, where p = (p x ,p )n \)' and q = {q x ,q y ,\)', then //// = 
</, where //=(/*/*,/* py, 1)' and qf = ( / * ?x, / * fly, 1)'. Thus, if B k is a 

10 transformation between two low resolution images I k and J h then B k is an 
approximation of the transformation between images 1 and J. That is, 5jf ~ 5 and dB k 
« J2?//r, where denotes the displacement values 

dB(x^) = (dx*dy 9 \) = {x-u*y-v 9 1), 5(xjU) = (w,v,l) (14) 
The quality of the approximation depends on the size of the factor k because of the 

1 5 inevitable loss of information in low resolution images. 

The multiresolution registration process starts by selecting a scale factor k 0 > 
1, where k$ is not too large, so that the scaled down images contain enough 
information to calculate a meaningful, although coarse, approximation of the 
transformation between the two images. Denote by / and J the target and reference 

20 full scale (input) images, respectively. Let T ku be the transformation between I ko and 

J k . Tf* is the first coarse approximation for the transformation T between / and J. 
Denote I k " = if* (I) . /*" is the result of applying the first approximation to the input 
target image. Now the approximation is refined by applying the method to and 
J k , where *, < k {h to yield the transformation T i% . The transformations 7^" and T k \ x 
25 are composed to yield a refined affine mapping M kl between the original images. 
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This process continues iteratively, such that at each step T k is calculated and used to 
refine M kj : 

M k ' = (15) 

Note that T k is calculated between I k k '" x and J k , where: 

5 A=M A '(/) (16) 

Because the affine transformations are closed under composition, M kj is an affine 

transformation* and is a direct mapping of the input image, so that the iterative 
process does not accumulate errors. The reduction factors k t can be selected 
arbitrarily. Preferably , however, the binary reduction, k, = k iA l2, />0, is used. This 

10 continues until the original resolution is restored. Figure 1 is a diagram of the data 
flow of multiresolution registration. 

The multiresolution method improves the performance of the registration in 
terms of the initial misalignment of the source and target images. Figures 2 A through 
2F present an example of the use of multiresolution to improve alignment. Figure 2A 

15 is a reference image. Figure 2B is a target image to be warped to match the reference 
image of Figure 2A. Note that Figure 2B has a small square in its upper left hand 
corner for reference. Figure 2C shows the results of warping Figure 2B using the 
gradient method without multiresolution. The normalized cross-correlation between 
Figures 2A and 2C is only 0.82. Figure 2D shows the results of warping Figure 2B 

20 using multiresolution. The normalized cross-correlation between Figures 2A and 2D 
is 0.93. This multiresolution registration was done in three steps. Figure 2E shows 
the image of Figure 2B after one multiresolution step, with a normalized cross- 
correlation of 0.8 1 with Figure 2A. Figure 2F shows the image of Figure 2B after two 
multiresolution steps, with a normalized cross-correlation of 0.91 . - 

25 Culling 

The registration method described above is global in the sense that it considers 
all the pixels of the source and target images. The solution of the equation system is 
based on the least squares criterion, in which all the pixel values contribute to the 
solution. Some of the pixels are irrelevant to the solution, as the corresponding pixels 

30 in the other image bear no resemblance to them, or may not even exist. The "vote" of 
these pixels not only cannot contribute, but may even add noise to the equation 
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system.. The present invention includes a rejection process that culls away pixels that 
are estimated to not contribute to the correct solution. Specifically, a fast local 
similarity measure is used: 

E{x,y) = X.Z,|U W + *> v + 30- Is ~ A (" + ^ v + y) + A| 0 7 ) 

5 where 7 V and /, are the average values in some local neighborhood around I s (xj>) and 
7,(jc,y), respectively. A pixel at (xj/) whose similarity value E(x,y) is smaller than a 
given threshold is rejected. . Determining an appropriate local neighborhood and an 
appropriate threshold is straightforward -for one ordinarily skilled in the art. 

Localization 

10 In most real life cases it is quite rare to find a global polynomial that 

transforms one image to another. In most cases the transformation characteristic is 
more general and nonuniform, with different variations over the image. It therefore is 
necessary to construct the transformation by an aggregation of many local 
transformations, which define a general elastic transformation. The present invention 

15 finds local transformations using the above global registration method over small 
windows and then interpolates the local transformations to define the global 
transformation. The localization permits using transformations of low degree, whose 
combination can define a global elastic transformation of high degree. The present 
invention preferably uses local affine transformations, taking advantage of the fact 

20 that affine transformations are closed under composition. 

However, there is a problem of working directly on small areas. A small 
window in one image is not the transformation of the corresponding window in the 
other image. This means that it is necessary to deal with partial matching. Moreover, 
it may occur that images of some windows are mismatched while the entire image 

25 matches. Thus, an independent registration between the windows most likely will 
fail. This suggests the necessity of defining a progressive process in which the 
registration is defined by a series of refined transformations. In the following, an 
hierarchical localization and a (weighted) interpolation are introduced. The process is 
illustrated by the example of Figures 3A through 3E. 

30 The localization process starts with the definition of a global transformation 

between the initial two images, the reference image of Figure 3 A and the target image 
of Figure 3B. The target image of Figure 3B is mapped to an intermediate image, 
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Figure 3C, by this global transformation . The intermediate image of Figure 3C and 
the reference image of Figure 3A are each partitioned into four windows. Then the 
corresponding windows are registered. This process continues recursively in each 
window 7 and is iterated down to a predefined level of localization. 
5 At localization level k, the images are partitioned into 2 A by 2* windows. Let 

L) . denote the local transformation centered at the (/,/)-window of the A-th level. 
First, the global transformation 7^ is computed between the input images I and J. The 
transformed image / = 7°(7) and the referee image J each are partitioned into four 
images: 1^, /<'„, /,' 0 , /,', and Jdo> ^» J lo> ^i'm respectively. The four local 
10 transformations Zj . are computed between the corresponding four pairs of images. 
The local transformations are then composed with the global transformation to yield 
T.y , the local composed transformations: 

A global transformation M x is defined as the weighted average of these four 
15 transformations: 

a/ 1 =y l y l w l .y d9) 

where the w tJ - are weights defined below. 

Figure 4 shows the data flow of the localization process. In general, at level /c, 

each local transformation L h } ) is composed of the upper level transformation 

(20) 

and the global transformation il/ is defined by 

M" =Y*~ x Y?~'y> u T (21) 

The transformed image / = and the reference image J are then partitioned again 
and the process is repeated. 
25 The local composed transformations are combined to form a global elastic 

mapping by an inverse distance weighted interpolation: 

where satisfies: 
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3. m-(x.^)>0 

and (x, v j/ /V ) is the sub-image center. Shepard (D. Shepard (1968), A two-dimensional 
5 interpolation function for irregularly spaced data, Proc. 23rd Natl. Conf. Of the ACM, 
ACM Press, New York, pp. 517-524) suggested the following function: 

, , d{{xMx LnyiJ ))- 

for any (xy) * (x />y j/ /f/ ) and 1 otherwise. The Euclidean distance between points a and 
b is denoted by d{a,b), and the value of the exponent u determines the smoothness of 
10 the transformation. Preferably, u=2 to emphasize the local influence. 

Note that the global transformation A/ is neither an affine transformation nor a 
linear transformation, but a general elastic transformation. The computation of each 
local transformation L k u is computed with the multiresolution and culling techniques. 
The example in Figures 3 A through 3E shows the importance of localization in 
15 cases where it is necessary to register images between which no affine transformation 
exists., for example, in cases where the deformation is local to a small part of one of 
the images. The global warped image of Figure 3C has a normalized cross-correlation 
with the image of Figure 3A of only 0.17. Figure 3D shows the results of using four 
local registrations (one level of localization): the normalized cross-correlation 
20 between the images of Figures 3A and 3D is 0.64. Figure 3E shows the results of 
using sixteen local registrations (two levels of localization): the normalized cross- 
correlation between the images of Figures 3 A and 3E is 0.81 . 

Although the preferred mode of localization is by halving the linear dimension 
of the subwindow at every localization level, as described above, the scope of the 
25 present invention includes all ways of subdividing the reference and target images into 
nonoveiiapping subimages, for example, dividing the linear dimensions of the 
subwindow by an integer greater than 2 at every localization level. 

Extension to Higher Dimensions 
As was mentioned above, the sampled ultrasonic volumes are deformed during 
30 acquisition. Because unavoidable deformation is inherently a 3D deformation, the 
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warp transformation that can align the deformed volume must be a true 3D warp. 
Otherwise, it would be like trying to register two slices lying on different planes. 

The extension of the gradient method to 3D is straightforward. All that is 
required is to treat the third dimension analogously to the other two dimensions. 
5 Define a 3D affine transformation: 
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Assume, as in the 2D case, that / and ./ are smooth differentiate functions 
from R 3 to R. What is needed is an affine transformation T{x x y,z) - (s,t,w>) such that 

l{xy,z) = J(s 9 t,w). Define the displacements dx-s - x s dy-t- and dz = w -z, so: 
1 0 Kxyj) = Ax + dx,y + dy\ z + dz) (25) 

If the displacement vector (dx,dy,dz) is small enough to be approximated by 
the linear term of a Taylor expansion, then 

J{x + dx,y + dya + dz) = J(x, v,z) + dxd v J{x,y y z) + cfyd v J(x,y,z) + dzdj(x,y,z) 

(26) 

15 Because 

dx = u - x ~ a [{ + {a X2 - \)x + a^y + a M z 
*/y = v - j; = a 2 i + tf 22 x + ("23 - 1 )j ; + ^ 

cfe = w - z = fl 3 1 + a 22 x + + (a 34 - 1 )z (27) 
it follows that 

20 / (x,y,z) = J(x 9 y, z) + («,,+ (« t2 - l)x + a l3 j; + a ]A z)d x J(x,y,z) 

+(a 2l + 22 x + (a 23 - 1);> + a 24 z)a i) J(x, < y ? z) + (a 3l + « 32 x + a„ v + («„ - l)z)a_ J(x, 

(28) 

Because it is assumed that Equation (28) holds for most of the voxels, there are 
obtained, as in the 2D case, a large number of equations which are solved by the 
25 pseudo inverse method, similar to the 2D case (Equation (9)). The multiresolution, 
culling, and localization techniques are extended and applied in a similar way. In 
particular, localization is done using 3D nonoverlapping subimages, preferably by 
partitioning the images initially into eight subimages, and, in each succeeding level of 
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localization, partitioning each of the subimages of the previous level into eight 
subimages. 

This formalism is easily extended by one skilled in the art to images of 
dimensionality higher than 3, and the scope of the present invention includes the 
5 registration of images of dimensionality higher than 3. In the localization step for 
images of dimensionality D, if the subimages are obtained by dividing the linear 
dimensions of the parent images by an integer N at every localization level, then each 

parent image is divided into N 2 ' subimages at every localization level. 



The gradient method defines an elastic transformation which warps the target 
volume to agree with the geometry of the reference volume. The method relies on the 
assumption that the two volumes are relatively close. Given two volumes, it is 
necessary to recover their relative positions. With 2D images this can be done 

15 manually, but it is almost impossible to have the same success in 3D because of the 
lack of proper 3D user interfaces. 

A 3D rigid transformation is a function of six variables, three for translation 
and three for orientation. Unlike the elastic gradient method, here an indirect method, 
i.e., a search method, is used. Because the two volumes are of the same type, a cross- 

20 correlation can be used as a similarity measure. The search for a global maximum is 
carried out in a six dimensional domain 

The problem of finding a global maximum for a multidimensional function is 
well known. This is a classic problem of optimization in which the main difficulty is 
to avoid converging to some local maximum. There are many solutions for non linear 

25 systems, such as Newton or Taylor techniques, which assume an initial proximity to 
the global solution, because they do not distinguish between global and local extrema. 
The preferred method in the present invention is the well established technique of 
"simulating annealing", which can elude local extrema. The parameters of the 
technique require data-dependent tuning, but this can be done easily by one skilled in 

30 the art. A simulating annealing method can also be applied to search for a non-rigid 
transformation. However, because search methods are sensitive to the number of 
parameters, the search time is impractical. 



10 
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Seamless Mosaicing 
After the coarse rigid alignment has been applied, the location of the target 
volume in terms of the reference volume coordinate system is known, as shown in 
Figure 5. This defines a clipping box contained within the overlapping portion of the 
5 two volumes. The target volume is warped by the 3D gradient method based on the 
context of two subvolumes } one from each volume, bounded by the clipping box. 
These two steps, i.e., rigid alignment and the 3D elastic warp, are the first two stages 
along the mosaicing pipeline illustrated in * in ire 6. Equation (22) is assumed valid 
throughout the target volume, not just within the clipping box. The rest of the 

10 pipeline stitches the warped volume over the reference volume. 

The warped target volume has to be composed over the reference volume in a 
seamless way. In most cases, adjacent volumes have different intensities, mainly 
because of directional variations. Some parts in one volume bear no resemblance to 
the corresponding part of the other volume or often do not exist. This requires the 

15 application of a 2D registration of the two slices at the seam face to further refine the 
misalignments at the connecting planes, and to blend (cross-dissolve) the two volumes 
to blur away the seam faces. 

The determination of the seam planes is not as simple as it may seem. The 
seam plane should pierce the clipping volume but not necessarily in the middle. 

20 Because we apply a 2D registration along this plane, it should contain significant 
information for a valid registration. A simple heuristic may be applied to look for the 
slice parallel to the volume face which has the most meaningful pixels, »where a 
meaningful pixel is defined as a pixel whose close neighborhood has a variation above 
some predefined threshold. Note, that there may be more than one seam plane. In 

25 many cases there are three seam planes, and if the target volume is set deep into the 
reference volume there may be up to six seam planes. 

Once the seam planes are determined, two slices are extracted, one from each 
volume for each seam plane. The images of the slices are registered by the 2D version 
of the elastic multiresolution gradient method described above. The 2D registration 

30 defines a set of transformations Mj (Equation (21)) one for each scam plane. Then, the 
target volume is deformed by a weighted blend of the 2D transformations defined as 
follows. Let p be an arbitrary point in the target volume and;?, the closest point (the 
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orthogonal projection) on the /-th seam plane. Let iv, be an inverse power of the 
Euclidean distance between p and p t : 



where the exponent u controls the distance decay. Then 



The blended registration function M improves the geometric registration in the 
proximity of the seam planes and reduces the distortions away from them. The gray 
levels of the two images are blended by a similar weighted sum to give the final gray 
level at point p. 

While the invention has been described with respect to a limited number of 
embodiments, it will be appreciated that many variations, modifications and other 
applications of the invention may be made. 




(31) 
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WHAT IS CLAIMED IS: 

1. A method for warping an input target image to match an input 
reference image, both images being at an original level of resolution, both images 
being composed of a plurality of pixels, the method comprising the steps of: 

(a) convolving the input target image with a low pass filter, thereby 
obtaining a filtered target image; 

(b) convolving the input reference image with said low pass filter, thereby 
obtaining a filtered reference image; 

(c) subsampling said filtered target image, thereby obtaining a subsampled 
target image; 

(d) subsampling said filtered reference image, thereby obtaining a 
subsampled reference image; 

(e) determining an affine transformation that maximizes a normalized 
cross-correlation between said subsampled target image and said 
subsampled reference image; and 

(f) scaling said affine transformation to apply to the original level of 
resolution, thereby obtaining a scaled factor of said affine 
transformation. 

2. The method of claim I, wherein said low pass filter is a box filter. 

3. The method of claim 1, wherein said affine transformation has a 
plurality of coefficients, and wherein said determination of said affine transformation 
includes a least squares solution for said coefficients. 

4. The method of claim 1, wherein said convolving, subsampling, and 
affine transformation determination are iterated at a plurality of levels of increasing 
resolution, wherein said subsampled target image, obtained at one of said levels, is 
used as the target image in a next of said levels, the method further comprising the 
step of: 
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(g) composing said scaled factors to obtain a composed affme 
transformation. 



5. The method of claim 1, further comprising the steps of: 

(g) computing a similarity measure between at least some the pixels of the 
input target image and at least some of the pixels of the input reference 
image; and 

(h) culling pixels for whicn said similarity measure is less than a threshold. 



6. A method for warping an input target image to match an input 
reference image, both images being of a dimensionality D, both images being at an 
original level of resolution, both images being composed of a plurality of pixels, the 
method comprising the steps of: 

(a) subdividing the input target image into a plurality of target subimages; 

(b) subdividing the input reference image into a plurality of reference 
subimages, there being a one-to-one correspondence between said 
reference subimages and said target subimages; 

(c) for each of said target subimages, determining a local transformation 
that warps said target subimage to match said corresponding reference 
subimage. 

7. The method of claim 6, wherein said determining of said local 
transformation is effected by the steps of: 

(i) convolving said target subimage with a low pass filter, thereby 
obtaining a filtered target subimage; 

(ii) convolving said reference subimage with said low pass filter, thereby 
obtaining a filtered reference subimage; 

(iii) subsampling said filtered target subimage, thereby obtaining a 
subsampled target subimage; 

(iv) subsampling said filtered reference image, thereby obtaining a 
subsampled reference subimage; 
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(v) determining a local affine transformation that maximizes a normalized 
correlation between said subsampled target subimage and said 
subsampled reference subimage; and 

(vi) scaling said local affine transformation to apply to the original level of 
resolution, thereby obtaining a scaled factor of said local affine 
transformation. 

8. The method of claim 7, wherein saiu iocal transformation is said scaled 
factor of said local affine transformation. 

9. The method of claim 7, wherein said convolving, subsampling, and 
affine transformation determination are iterated at a plurality of levels of increasing 
resolutions, wherein said subsampled target subimage, obtained at one of said levels 
of resolution, is used as said target subimage in a next of said levels of resolution, the 
method further comprising the step of: 

(vii) composing said scaled factors to obtain said local transformation. 

10. The method of claim 7, wherein said determining of said local 
transformation further includes the step of: 

(vii) computing a similarity measure between at least some the pixels of 
said target subimage and at least some of the pixels of said 
corresponding reference subimage; and 

(viii) culling pixels for which said similarity measure is less than a threshold 

11. The method of claim 6, wherein both said subdividing, and said 
determining of said local transformation, are performed recursively at a plurality of 
levels of localization, the method further comprising the steps of: 

(d) for each of said target subimages, composing said local transformation 
with a prior transformation, thereby obtaining a local composed 
transformation; 

(e) forming a weighted sum of said local composed transformations, 
thereby obtaining a global transformation; and 
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(f) applying said global transformation to said target image, thereby 
obtaining a transformed target image; 
said transformed target image, obtained at one of said levels of localization, being 
used as the input target image in a next of said levels of localization. 

12. The method of claim 1 1 . wherein, in each of said levels of localization 
after a first of said levels of localization, said prior transformation is said local 
transformation from a preceding level of localization. 

13. The method of claim 6, wherein the input target image is divided into 

;Y 2 target subimages and wherein the input reference image is divided into N~ 
reference subimages, where N is an integer greater than 1 . 

14. The method of claim 1 3 , wherein N is 2. 

15. The method of claim 6, further comprising the step of: 

(d) rigidly transforming the input target image to match the input reference 

image. 

16. The method of claim 15, wherein said rigid transformation is effected 
by steps including maximizing a normalized cross-correlation between the input target 
image and the input reference image. 

17. The method of claim 16, wherein said maximizing is effected by 
simulated annealing. 

18. A method for blending a three-dimensional target image with a three- 
dimensional reference image, the target image including a plurality of points, the 
method comprising the steps of: 

(a) determining at least one seam plane shared by the target image and the 
reference image; 

(b) for each of said at least one seam plane: 
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(i) extracting an input target slice from the target image along said 
seam plane, 

(ii) extracting an input reference slice from the reference image 
along said seam plane, 

(iii) subdividing the input target slice into a plurality of target 
subslices, 

(iv) subdividing the input reference slice into a plurality of 
reference subslices, there being a one-to-one correspondence 
between said reference subslices and said target subslices, 

(v) for each of said target subslices, determining a local 
transformation that warps said target subimage to match said 
corresponding reference subimage, 

(vi) recursively performing both said subdividing and said 
determining of said local transformation at a plurality of levels 
of localization, said recursion further including the steps of: 

(A) for each of said target subslices, composing said local 
transformation with a prior transformation, thereby 
obtaining a local composed transformation, and 

(B) forming a weighted sum of said local composed 
transformations, thereby obtaining a global 
transformation; 

(c) applying a weighted sum of said at least one global transformation to at 
least one of the points. 
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traasformation and a fine elastic warp. The fine transformation assumes that the two volumes are close enough for the mapping between 
them to be expressed in terms of their gradient values. 



FOR THE PURPOSES OF INFORMATION ONLY 



Codes used to identify States party to thePCT on the front pages of pamphlets publishing international applications under the PCT. 



AL 


Albania 


ICS 


Spain 


LS 


L^sotlto 


SI 


Slovenia 


AM 


Armenia 


FI 


Finland 


LT 


Lithuania 


SK 


Stovakia 


AT 


Austria 


VR 


France 


LU 


Luxembourg 


SN 


Senegal 


AU 


Australia 


GA 


Gabon 


LV 


Latvia 


sz 


Swaziland 


AZ 


Azerbaijan 


Gn 


United Kingdom 


MC 


Monaco 


TD 


Chad 


BA 


Bosnia and Herzegovina 


GE 


Georgia 


MD 


Republic of Moldova 


TG 


Togo 


BB 


Barbados 


Gil 


Ghana 


MG 


Madagascar 


TJ 


Tajikistan 


BE 


Belgium 


GN 


Guinea 


MK 


The former Yugoslav 


TM 


Turkmenistan 


BF 


Burkina Faso 


GR 


Greece 




Republic of Macedonia 


TR 


Turkey 


BG 


Bulgaria 


1IU 


Hungary 


ML 


Mali 


'IT 


Trinidad and Tobago 


BJ 


Benin 


)E 


Ireland 


MN 


Mongolia 


UA 


Ukraine 


BK 


Brazil 


IL 


Israel 


MR 


Mauritania 


UG 


Uganda 


BY 


Belarus 


IS 


Iceland 


MW 


Malawi 


US 


United States of America 


CA 


Canada 


IT 


Italy 


MX 


Mexico 


uz 


Uzbekistan 


CF 


Central African Republic 


JP 


Japan 


NE 


Niger 


VN 


Viet Nam 


CG 


Congo 


KE 


Kenya 


NL 


Netherlands 


YU 


Yugoslavia 


CH 


Switzerland 


KG 


Kyrgyzslan 


NO 


Norway 


ZW 


Zimbabwe 


CI 


Cote d'lvoire 


KP 


Democratic People's 


NZ 


New Zealand 






CM 


Cameroon 




Republic of Korea 


PI. 


Poland 






CN 


China 


KK 


Republic of Korea 


FT 


Portugal 






CU 


Cuba 


KZ 


Kazakstan 


RO 


Romania 






C2 


Czech Republic 


LC 


Saint Lucia 


RU 


Russian Federation 






DE 


Germany 


LI 


Liechtenstein 


SD 


Sudan 






DK 


Denmark 


LK 


Sri Lanka 


SE 


Sweden 






EE 


Estonia 


LR 


Liberia 


SG 


Singapore 







INTERNATIONAL SEARCH REPORT 


International application No. 




PCTAJS97/20950 



A. CLASSIFICATION OF SUBJECT MATTER 
1PC(6) : G06TC 9/62 

USCL : 382/215 " . ^ . Jin n 

According to International Patent Classification (IPC) or to both national classification and IPC 



a FIELDS SEARCHED 

Minimum documentation searched (classification system followed by classification symbols) 

U.S. : 382/215, 294, 293, 291, 284, 216, 209; 358/450; 395/136 
Documentation searched other than minimum documentation to the extent that such documents are included tn the fields scarche* 

Electron r^i* base consulted during the international search (name of data base and, where practicable, search terms used) 
IPSS database search terms: transform, nffine, blend, match, register, filter, convolve, subsample, seam 



C DOCUMENTS CONSIDERED TO BE RELEVANT 



Category* 



A,P 



Citation of document, with indication, where appropriate, of the relevant passages 



US 5,657,402 A (BENDER ET AL) 12 AUGUST 1997, ENTIRE 
DOCUMENT 

US 5,633,951 A (MOSHFEGHI) 27 MAY 1997, ENTIRE 
DOCUMENT 

US 5,611,033 A (P1TTELOUD ET AL) 11 MARCH 1997, ENTIRE 
DOCUMENT 

US 5,550,937 A (BELL ET AL) 27 AUGUST 1996, ENTIRE 
DOCUMENT 



Relevant to claim No. 



1-18 
18 

1-18 
1-18 



| [ Further documents are listed in the continuation of Box C. Q See patent family annex. 



Special categories of cited documentor 

document defining tba general iUU of the art which u not oonsioWi 
to bt of particular reference 

earlier docuaent published on or after Um international filing data 

document which may throw doubts on priority elaim(s) or which ia 
ciud to establish the publication date of another citation or other 
special reaaon (aa ■pacified) 

document referring to an oral dbcloaura. um. exhibition or other 



later document published after the international filing date or priority 
date and not in conflict with the application but cited to understand 
the principle or theory underlying the invention 

document of particular relevance; the claimed invention cannot be 
considered novel or cannot be considered to involve an inventive step 
when the document is taken alone 

document of particular relevance; the claimed invention cannot be 
conaidared to involve an inventive step when the document is 
combined with one or more other auch documents, such combination 
being obvious to a person skilled in the art 



■p* document published prior to the mUmational filing date but later than document member of the same patent family 


Date of the actual completion of the international search 
13 FEBRUARY 1998 


Date of mailing of the international search report 

j^l 9 /; BAY 1998 


Name and mailing address of the ISA/US 
Commissioner of Patents and Trademarks 

Box PCT 

Washington, D.C. 20231 
Facsimile No. (7<> 3 ) 305-3230 


AuthortzfisTofpgc^^ 

MICHAEL RAZAV1 
Telephone No. (703) 305^713 



Form PCT/ISA/210 (second sheetX-July 1992)* 



